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ABSTRACT 

We describe a phenomenological model for molecular hydrogen formation suited for applications in 
galaxy formation simulations, which includes non-equilibrium formation of H2 on dust and approxi- 
mate treatment of both its self-shielding and shielding by dust from the dissociating UV radiation. 
The model is applicable in simulations in which individual star forming regions - the giant molecular 
complexes - can be identified (resolution of tens of pc) and their mean internal density estimated 
reliably, even if internal structure is not resolved. In agreement with previous studies, calculations 
based on our model show that the transition from atomic to fully molecular phase depends primarily 
on the metallicity, which we assume is directly related to the dust abundance, and dumpiness of 
the interstellar medium. The dumpiness simply boosts the formation rate of molecular hydrogen, 
while dust serves both as a catalyst of H2 formation and as an additional shielding from dissociating 
UV radiation. The upshot is that it is difficult to form fully-shielded giant molecular clouds while 
gas metallicity is low. However, once the gas is enriched to Z ~ 0.01 — 0.1Z©, the subsequent star 
formation and enrichment can proceed at a much faster rate. This may keep star formation efficiency 
in the low-mass, low-metallicity progenitors of galaxies very low for a certain period of time with the 
effect similar to a strong "feedback" mechanism. The effect may help explain the steep increase of 
the mass-to-light ratio towards smaller masses observed in the local galaxy population. We apply the 
model and star formation recipes based on the local amount of molecular gas to an output from a 
cosmological simulation of galaxy formation and show that resulting global correlations between star 
formation and gas and H2 surface densities are in good agreement with observations. 
Subject headings: cosmology: theory - galaxies: evolution - galaxies: formation - stars:formation - 
methods: numerical 



1. INTRODUCTION 

Detailed understanding of galaxy formation remains 
one of the most important goals of modern astrophysics. 
A critical challenge for achieving this goal, in addition 
to a sheer variety of physical processes and the range of 
spatial and temporal scales involved, is our poor under- 
standing of how gas is converted into stars under differ- 
ent conditions. This conversion controls critical aspects 
of the final galaxy properties from its total stellar mass 
and star formation history to its morphology, which de- 
pends on the relative timing of gas conversion to stars 
and the epoch of major mergers and rapid mass assem- 

Recent observational and theoretical progress in un- 
derstanding star formation on the scale of molecular 
clouds points towards a prescription for star formation 
on larger, galactic scales. In particular, the observational 
data support the theoretically motivated view of molec- 
ular clouds in which the specific star formation rate per 
local free fall time in the molecular gas is essentially in- 
dependent of the local gas den s ity dKrumholz fc McKe 
l2005t iKrumholz fc^ranM2007at iKrumholz et alJ 12008 " 
This is fortunate because the efficiency of star formation 
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can then be modeled in a simulation as long as (i) one 
can resolve and identify the regions corresponding to star 
forming molecular clouds in a simulation and (ii) achieve 
resolution sufficient to model their mean internal density 
properly. 

Until recently, these conditions were not easily 
achievable in cosmological simulations due to numerical 
limitations on spatial and mass resolution. The stan- 
dard approach in galaxy formati on mo deling so far (c.f. 
recen t studies bv iBrooks et al.1 120071: iGovernato et al. | 
2007 : ICeverino fe Klvpinl 20071: IScannapieco et al.l 
20081 ICroft et all 120081: lOppenheimer fe Pavel 120081: 
Tasker fe Brvanl l2008t iMaver et all I2008L and many 



earlier works, comprehensively reviewed in the last 
reference) was to adopt a recipe which ties the local 
star formation rate density to the gas density via a 
universal relation, often with threshold conditions on 
the gas properties (density, temperature, etc). Such 
recipes are loosely based on the empirical correlations 
observed for local gal axies (the Kennicu t t-Schmidt, 
herea f ter KS , relations: ISchmidtJll959l Il963t iKennicuttl 
119891 Il998bf ). However, these relations have only 
been studied relatively well for nearby massive or star 
bursting galaxies. There is, however, a growing body of 
evidence that the relations for galaxies of lower surface 
brightness and/or metallicity may be quite different 
(see, e.g., iRobertson fc Kravtsovl [20081 fo r references 
and discussion and the recent study by iBigiel etlrtl 
12001 . 

Ind ee d, as we emph a size b el ow (see also lElmegreenl 
19891 [1991: iSchavd 120011 iPelupessv et al.l 120061 : 
Krumholz et al. I I2008bf ). the transition from atomic 
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to fully molecular gas depends strongly on the local 
gas metallicity and dissociating UV flux. Thus, for 
example, we expect that in low metallicity galaxies gas 
becomes fully molecular (and thereby conducive to star 
formation) at higher gas density compared to more 
massive, higher metallicity systems. This means that on 
the global scale of a galaxy gas may be converted into 
stars at a slower rate per unit mass of gas or, in other 
words, the star formation efficiency is lower, even if the 
density distribution of ISM is the same. Not only is 
this regime applicable to dwarf or low surface brightness 
galaxies today, but also to the majority of progenitors 
of massive galaxies at high redshifts. 

The lower efficiency of gas conversion into stars at high 
redshifts may have profound implications for galaxy evo- 
lution from shaping the faint end of the luminosity func- 
tion to the morphological mix of galaxies. Tradition- 
ally, in cosmological simulations and semi-analytic mod- 
els the efficiency of star formation in < L» galaxies is sup- 
pressed by supernova feedback. However, the efficiency 
of such suppression is not clear and in many observed 
dwarf galaxies star formation is inefficient without ob- 
vious signs of active feedback. We argue that a similar 
suppression can be provided by the inherent difficulty 
of building self-shielding molecular regions within low- 
density, low-metallicity ISM of smaller galaxies. 

Current self-consistent cosmological simulations 
of galaxy formation can resolve s cales of molec- 
ular clouds (i.e., t ens o f parsecs, IKravtsovl 20031 : 
Ceverino fe Klvpinl [20071: iTasker fc Brvanl 120081 : 



Razoumovl l2008t iGibson et a l. 2008) and can, therefore, 
explore the effect of ISM process on the efficiency of 
star formation self-consistently. Observational evidence 
for the universal efficiency of gas conversion per free 
fall time in molecular clouds can then be adopted as a 
well-motivated and robust star formation prescription, 
even if the internal density structure of the star forming 
regions is not well resolved. Such recipe is promising 
but requires a model for the formation and destruction 
of molecular gas in the simulations in order to correctly 
identify the regions of current star formation. 

Several such models have recently bee n developed and 
implemented in the SPH simulations dPelupessy et al.l 
20061 iBooth et alj [20071: iRobertson fc Kravtsovl l200g T 
Pelupessv et al.l ( 20061 ) have followed reactions of H2 for- 
mation and destruction, following reactions of H2 for- 
mation on dust grains and dust and self-shielding from 
UV radiation and using a sub-grid model for gas clouds, 
which uses observed cloud scaling relations. One of the 
main conclusions of their study was the strong depen- 
dence of the molecular gas fracti on on the ambient metal- 
licity of ISM. lBooth et all (|2007t ) have adopted a different 
approach, in which they model formation of molecular 
clouds using a model based on collisions between clouds 
representing molecular clouds as ballistic particles coagu- 
lating upon collision. They have demonstrated that with 
such a model many observed p roperties of disc galax- 
ies can be reproduced. Recently, IRobertson fc Kravtsovl 
(2008) presented a model of star formation based on 
the local molecular gas content. The latter was mod- 
eled by using pre-computed tables of molecular fraction 
as a function of local gas properties (density, temper- 
ature, metallicity and UV flux). The authors showed 
that the model implied a significant steepening of the 



Kennicutt-Schmidt relations in low surface density, low 
molecular fraction environments characteristic of dwarf 
and low surface brightness galaxies, as well as of the out- 
skirts of disks in larger galaxies. This study has indi- 
cated that one of the key factors controlling the amount 
of diffuse molecular gas in the lower density ISM is UV 
radiation from the recent local star formation. 

In this paper we present a model of molecular gas evo- 
lution and Hi2-based star formation together with their 
implementation in the Adaptive Refinement Tree (ART) 
code for cosmological simulations. The ART code is 
a Eulcrian hydrodynamics+V-body code, which uses 
Adaptive Mesh Refinement (AMR) technique to reach 
high resolution in the regions of interest (i.e., in high- 
density regions of ISM in the case of this study). We 
present details of the model of molecular hydrogen forma- 
tion and disruption and test the model against available 
observations. We discuss several star formation recipes 
for the gas conversion in molecular clouds and present 
resulting Kennicutt relations, when such recipes and the 
molecular hydrogen model are applied to a realistic cos- 
mological simulation of a Milky Way progenitor. In a 
follow-up study, we plan to explore the differences of such 
model compared to the old star formation prescriptions 
in full cosmological simulations. 

2. METHOD 

2.1. Description of Simulations 

The cosmological simulation used as a testing ground 
for our model for molecu lar hydrogen form ation has 
been described in detail in iTassis et alj (|2008l ) as a run 
"FNEC-RT" . Here we only recap that the simulation was 
performed with the Eulerian, gas dyna mics + jV-body 
Adaptive Refinement Tree (ART) code (|Kravtsov et all 
Il997t IKravtsovl [19991: iKravtsov et~aTll2002D . which uses 
adaptive mesh refinement (AMR) in both the gas dynam- 
ics and gravity calculations to achieve a needed dynamic 
range. 

The simulation follows a Lagrangian region corre- 
sponding to five virial radii of a system, which evolves 
to approximately Milky Way mass (M rs 10 12 M Q ) at 
z = 0. The mass resolution in dark matter of mj m = 
1.3 x 10 6 M Q (and corresponding resolution in gas dy- 
namics of m gas = 2.2 x 10 5 M Q ). This Lagrangian region 
is embedded into a cubic volume of 6ft. -1 comoving Mpc 
on a side, which is resolved to a lesser extent throughout 
the simulation. 

The simulation includes star formation and super- 
nova enrichment and energy feedback, and uses self- 
consistent 3-D radiative transfer of UV radiation from 
indivi dual stellar particles using the OTVET approxima- 
tion (|Gnedin fc AbelfeOOlh . The simulation follows non- 
equilibrium chemical network of hydrogen and helium 
and non-equilibrium cooling and heating rates, which 
make use of the local abundance of atomic, molecu- 
lar, and ionic species and UV intensity, followed self- 
consistently during the course of the simulation. We use 
the z — 4 output of the simulation as at this epoch the 
gaseous disk of the main galaxy is relatively quiescent 
and has not experienced major mergers since z « 6. 

A new ingredien t in the simu l ation, beyond what has 
been described in ITassis et al.l ([2008D . is an empirical 
model for formation and shielding of molecular hydro- 
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gen on the interstellar dust, which is described in the 
following section. 

2.2. The model for tracking molecular hydrogen 

The two key processes controlling the abundance of 
molecular gas in the ISM enriched with metals are the 
formation of H2 on dust grains and its shielding from 
photodissociation by the interstella r radiation field by 
itself (self-shielding) and by dust (see lGlover fc Mac Low! 
2007i for a recent review and references to a large body 
of prior work) . 

Unfortunately, the exact treatment of shielding cannot 
be achieved in modern cosmological simulations, as it re- 
quires three-dimensional radiative transfer, including ra- 
diative transfer in the Lyman- Werner bands. Therefore, 
we adopt a phenomenological treatment of molecular hy- 
drogen shielding, which is tuned to reproduce available 
observations (see below). 

In the optically thin regime (no shielding) and in the 
absence of dust, equations for the balance of neutral 
(atomic) hydrogen and molecular hydrogen can be writ- 
ten as 

Xr i = R(T)n e Xu n — -^h iTh i 2Xh 2 , 

*H 2 =*g, (1) 

where Xi denotes a number fraction of baryons in species 

i, 

Xi = rii/nb, 

and Ub is the total number density of baryons. In equa- 
tions ([T]) R(T) is the recombination rate, n e is the num- 
ber density of free electrons, Thi is the hydrogen ioniza- 
tion rate, and the term X^ includes multiple processes 
of formation and destruction of molecular hydrogen in 
the gas phase reactions, via H~ and H2 + ions. The lat- 
ter reactions are compr e hensively summarized elsewhere 
(iShapiro fc Kand fl987t [Abel et all [19971: iGalli fc Pallal 
1998: iGlover fc JappsenlkoOrt IG1 over fe Abelll2008T and 
are not important for our purposes; we keep them in the 
equations to ensure that our model also works in the 
metal-free regime. 

Physically, the effect of the geometric shielding by dust 
grains and H 2 self-shielding can be described by effective 
shielding factors, Sd and Sii 2 , which take values from 
to 1. Incorporating these factors, and adding molecular 
hydrogen formation on dust, equations |T]) become 

Xm = R(T)n e Xau — S^XhiLhi — 2Xn 2 , 

Xn 2 = SdSn 2 Xft 2 + RdnbXni (Xni + 2Xh 2 ) , (2) 

where we assume that the dust absorption also affects 
atomic hydrogen ionization rate. In principle, this is not 
a crucial assumption since one would expect that hydro- 
gen is predominantly neutral in regions where dust ab- 
sorption is important. In practice, however, the OTVET 
approximation for modeling radiative transfer, which we 
use in our simulations, is somewhat diffusive. In sim- 
ulations presented here, this diffusion leads to ionizing 
radiation leaking into neutral, self-shielded gas. The nu- 
merical diffusion effect is not significant, but introducing 
a factor Sd in equation (J^l) allows to essentially elim- 
inate it altogether. In a completely analogous manner, 
we also reduce the photo-heating rate of the gas in the 
evolution equation for the gas internal energy. 



The specific form of the shielding factors Sd and Sh 2 
cannot be computed without a detailed treatment of 
radiative transfer in the Lyman- Werner bands. Our 
goal, therefore, is to find physically plausible functional 
forms that are motiv ated by the observ a tional data. To 
this end, we follow iDraine fc Bertoldil (|1996l see also 
IGlover fe Mac Low! 12001 " and define these two factors 
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and 
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where x = Ns 2 /(5 x 10 14 cm 2 ) and aa, e s and ^h 2 are 
treated as adjustable parameters. We find that we obtain 
the best fit with the observational data (as we demon- 
strate below) for values of these parameters cx^eff = 
4 x 10 1 cm and wh 2 = 0.2, which are slightly d ifferent 
from those adopted by IGlover fc Mac LowT (|2007f) . 

The formation rate coefficient of molecular hydro- 
gen on dust, Rd, has been studied extensively in the 
past. While reviewing comprehensively this large body 
of work is beyond the scope of this paper, we would 
like to note that a the o retica lly cal culated formula from 
iHollenbach fc McKed (fl979h and iBurke fc Hollenbachl 
(119831) has been most c ommonly used in the past (see 
ICazaux fc Spaansl I2004L for a recent review) . In this 
work, following our main ideology of using the most em- 
pirically determined quantities, we adopt an obse r yation - 
ally determined value for Rd from lWolfire et ail {2008). 
We also assume that the dust-to-gas ratio scales linearly 
with gas metallicity Z (in solar units) and that the gas is 
clustered on scales unresolved in our simulations, so that 
the effective H2 formation rate is higher in proportion to 
the gas clumping factor, C p = (p 2 )/(p) 2 , 



R d = 3.5 x 10- 17 ZC p cm 3 /s. 



(3) 



The clumping factor C p can be viewed as a parameter 
of our model. However, there is also some theoretical 
arguments for its plausible range of values. Numerical 
simulations of turbulent molecular clouds typically find a 
lognormal density distribution inside the clouds with the 
width that depends on the average Mach number of the 
flows: of no » 1-2 (|Padoan et all 119971 : lOstriker etall 
I200H iMcKee fc Ostrikerl l200l . Since for a lognormal 
distribution, 
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most likely values for the clumping factor are in the range 
C p ~ 3 — 10. Empirical evidence also indicates that gas 
in molecular clouds is highly clumped with the ratio of 
typical molecular gas densi ty to the mean cloud de nsity 
of - 30 - 100 (see §3.1.2 in lMcKee fc Ostrikerfl2007L and 
references therein). 

Finally, in order to specify the column densities inside 
the molecular clouds, we use a Sobolev-like approxima- 
tion: 

Ni re ni L Sob , (4) 



where 
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Fig. 1. — Comparison of the total hydrogen column density 
from the Sobolev-like approximation and the true column density 
as integrated along random lines of sight through the simulation 
box. Gray points show individual lines of sight; black solid line 
show the average iVg b iev for a given N^mc, while black dashed 
lines give the rms scatter. The thin black dotted line is a diagonal 
of the plot. 




Fig. 2. — The intensity of the interstellar radiation field at lOOOA 
as a functio n of gas densi ty in the parent simulation FNEC-RT, in 
units of the Draine (1978) field (10 photons/ cm 2 / s/ sterrad/ eV). 
Since the parent simulation does not include dust shielding, radi- 
ation at 100CL4 in that simulation was computed in the optically 
thin regime. 



We have verified that equation (T4|) provides an essentially 
unbiased estimate for the true column density (obtained 
by integrating along random lines of sight) within the 
range 3 x 10 20 cm 2 < 7V H i + 2N^ 2 < 3 x 10 23 cm 2 with a 
scatter of about a factor of 2, as illustrated in Figure [T] 

In the gas dynamics solver and in the calculation of 
the cooling function, we use the thermodynamically cor- 
rect relations between the gas internal energy, pressure, 
and temperature, without assuming a polytropic approx- 
imation, which becomes invalid for mostly molecular gas 
(Turk, Gncdin, & Abel, in preparation). 

This model is implemented in the ART code. Because 
individual molecular clouds are not resolved in our sim- 
ulations, the specific numerical implementation suffers 
from numerical artifacts due to finite spatial resolution. 
We discuss these artifacts and our approach for their mit- 
igation in the Appendix. 



2.3. Star Formation 

We explore three simple star formation recipes, based 
on the local density of molecular hydrogen, in order to 
test the dependence of global KS relations on the un- 
derlying assumptions. We can write the star formation 
rate in a molecular cloud in terms of an efficiency fac- 
tor e, equivalent to the inverse of the (molecular) gas 
consumption timescale, as 

P* = epH 2 , (5) 

where is the star formation rate per unit volume, and 
pn 2 is the molecular hydrogen density. We can then ex- 
press e in terms of a dimensionless effective efficiency eg 
per local free-fall time of the gas, so that 
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where r s f = rg = y / 3Tr/32Gp for a uniform sphere. 
Since stars form out of regions of molecular complexes 
with sizes not necessarily resolved in our simulations, 
the choice of a relevant star formation timescale relies 
on assumption about the structure of the gas in sub-grid 
scales. We have tried two different assumptions about 
r s f. 

SF1 the constant time scale r s f = rg (100cm -3 ), assum- 
ing all H2 in a grid cell resides in molecular clouds 
of average density » 100 cm -3 ; 

SF2 the time scale, which depends on the 
cell's gas density with an upper limit 
T sf = min[r ff (100cm~ 3 ), Tff^c,,;;)]. 

SF3 the time scale, which simply depends on the density 
of the cell r s f = Tg(p . ce u) (this recipe is simila r to 
the one used by iRobertson fc Kravtsovl (|2008| ) in 
their models of isolated galaxies). 



We adopt e ff = 0.01 (|Krumholz k Tanll2007bh in all of 
the calculations presented below; we also use a molecu- 
lar fraction threshold for computational efficiency: star 
formation is allowed only in cells which have a molecular 
fraction higher than fn 2 =0.1. The latter criterion has 
a negligible impact on the total star formation rate in 
a galaxy, but is highly beneficial from a computational 
point of view, avoiding formation of a large number of 
very small stellar particles in low density regions. 

3. RESULTS 

In this section we present results illustrating the de- 
pendence of the relative atomic and molecular gas abun- 
dances on metallicity, gas dumpiness, and UV flux. We 
also explore the correspondent dependencies of the global 
Kennicutt-Schmidt relations with different local star for- 
mation prescriptions described in the previous section. 

Using t he z = 4 output of the parent FNEC-RT sim- 
ulation of iTassis et alj (|2008f ) (see § I2.1|) as initial condi- 
tion, we run a simulation with the model for molecular 
hydrogen and star formation for lOOMyrs. The base- 
line simulation FNEC-RT does not incorporate shield- 
ing of molecular hydrogen, and thus no molecular clouds 
form in that simulation. The chosen interval of time 
of 100 Myr is considerably longer than the char acteris- 
tic lifetime of molecular clouds (~ 10 7 yrs, iBlitz &: Shul 



TABLE 1 
Simulations 



Simulation 


Metallicity 


Clumping 


UV flux 


SF 




Z/Zq 


Factor C p 




Recipe 


A 


1 


10 


1 


SF2 


Bl 


0.3 


10 


1 


SF2 


B2 


0.1 


10 


1 


SF2 


CI 


1 


1 


1 


SF2 


C2 


1 


100 


1 


SF2 


D 


1 


10 


10 


SF2 


E 


1 


10 


1 


SF1 



119801 : iBlitz et al.lf2007t ) and is therefore sufficiently long 
for our purposes. We have also verified that our results 
have converged: various distribution presented below are 
virtually indistinguishable if we chose a 50Myr time in- 
terval instead. 

The baseline simulation FNEC-RT also includes effects 
of both energy and metal feedback to the interstellar 
medium, to keep the tests clean both forms of feedback 
are turned off and the interstellar medium metallicity 
is kept at a constant value after our test simulation is 
started. 

In our fiducial simulation (A) the metallicity is kept 
constant at solar value, we use a clumping factor C p — 
10, and the star formation recipe used is SF2. The UV 
flux is taken directly from the parent simulation FNEC- 
RT and is shown in Figure [5J Our model galaxy has 
a significantly higher star formation rate and is smaller 
at z — 4 galaxy than the Milky Way today, so the in- 
terstellar radiation field is substantially higher than the 
canonical Drainc (1978) value adopted for the Milky Way 
ISM at the present time, but is consistent with the esti- 
mates of the inters tellar UV field in high redshift galaxies 
(|Chen et al.ll2008lh 

In each of the simulations used in our parameter study, 
we vary a single parameter away from its fiducial simu- 
lation value. Table [1] describes the values of these free 
parameters in each of the simulations discussed in this 
section. 

A visual representation of our fiducial simulation A 
is shown in Figure [3J The molecular gas traces spiral 
arms well, which are less pronounced in the atomic gas. 
The molecular disk is both smaller and thinner than the 
atomic disk. However, since our spatial resolution is only 
50 pc, we do not resolve individual molecular clouds. In- 
stead, the orange colored surface in the image shows the 
boundary of the mostly molecular (/h 2 > 0.5) gas. 

3.1. Gas phases and atomic-to-molecular transition in 
the ISM of model galaxies 

Figure |4] shows the temperature in degrees Kelvin plot- 
ted against the total volume number density of neutral 
hydrogen for our fiducial simulation (A). Points in this 
and subsequent plots in this subsection show cells at the 
highest resolution level (level 9; the physical size of 52 pc 
at 2 = 4), which cover the large fraction of the disks of 
galaxies forming in the high-resolution Lagrangian region 
of the simulation. This plot clearly shows a well devel- 
oped multi-phase structure of the ISM in these disks. 
In particular, three different gas phases are evident: (i) 
the hot, ionized, low-density gas in the upper left part 
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Fig. 3. — A face- and edge-on views of the largest galaxy in 
our fiducial simulation A. Blue translucent color represents volume 
rendering of atomic hydrogen density, while orange color shows the 
location of the molecular gas (as the isosurface of /h 2 = 0.5 value). 
This figure is best viewed in color. 




10" 2 1 10 2 10 4 
n H (cm -3 ) 

Fig. 4. — Temperature plotted against the total gas number 
volume density for all maximum refinement level (level 9) cells at 
z = 4, for simulation A. 

of the diagram; (ii) the warm neutral medium around 
T ~ 10 4 K and n H ~ 0.1 - 10 cm" 3 and (iii) the cold 
neutral medium with T ~ 10 — 100 K at tin > 100 cm -3 . 
The transition from the warm neutral to cold neutral 
phase occurs over a narrow range of gas densities around 
a few tens cm~ 3 , in good agreement with the models of 
the Milky Way interstellar medium (whi c h has similar 
metallicity to our model) of IWolfire etaT] (|2003l ). 

Although not apparent in Figure HI the cold neu- 
tral medium undergoes a transition from atomic to 
fully molecular phase at the gas density of nn ~ 
100 cm~ 3 . Figure [5] shows the molecular fraction, Jh 2 = 
2?ih 2 / (nHi + 2?iH 2 ), as a function of the total neutral 
hydrogen volume number density. The three panels cor- 
respond to simulations with decreasing metallicity (left 
panel: simulation A, Z = Z & ; middle panel: simulation 
Bl, Z = 0.3Z G ; right panel: simulation B2, Z = 0.1Z o ). 
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Fig. 5. — Molecular hydrogen fraction in level 9 cells, plotted against the total neutral gas number volume density, 
corresponds to the density of binned points on the plane. Left panel: simulation A (solar metallicity); middle panel: 
(metallicity 0.3 solar); right panel: simulation B2 (metallicity 0.1 solar). 



The grayscale 
simulation Bl 



A sharp transition from fully atomic to fully molecular 
gas occurs at high gas densities, with the characteris- 
tic density of the transition increasing with decreasing 
metallicity. The gas density at which the molecular frac- 
tion is 50% scales with metallicity approximately as 



nt ~ 30 — 



cm 



(7) 



This trend is not surprising, as the atomic to molecu- 
lar transition is facilitated by the dust opacity: when 
dust opacity is sufficiently high to shield the gas from 
UV radiation, the transition to molecular medium can 
occur (self-shielding only becomes dominant at higher 
molecular fractions). As the metallicity increases, the 
amount of dust available to shield the gas from UV also 
increases, and the atomic-to-molecular transition can oc- 
cur at lower densities. 

To compare these results with observational measure- 
ments in Figure Owe plot the molecular fraction against 
the total neutral hydrogen column density, using a log- 
arithmic vertical scale, for simulations A, Bl, B2. We 
overplot the FUSE measurements of the molecular frac- 
tion in the MW halo and disk (left panel), in the LMC 
(m iddle panel), an d in th e S MC (right panel ) com piled 
bv lBrowning et al.l (|2003f ) and lGillmon et all (|2006l ). In 
the left panel we al so plot data for the Milky Way from 
IWolfire et ail (120081) , and on the ri ght p anel we show the 
data for SMC from iLerov et al.l (|2007f ). The metallic- 
ities of the ISM in these galaxies approximately corre- 
spond to the metallicities adopted in the models against 
which we compare them. In all cases, the column den- 
sity at which the molecular fractions begin to increase 
rapidly and the sharpness of the transition are well repro- 
duced in our simulations. For the solar metallicity case, 
this transition corresponds to the maximum surface den- 
sity of atomic hydrogen of about IOM0/PC 2 , in concor- 
dance with the measurements of molecul ar and atomic 
gas surface densities in nearby galaxies (Wong fe Blitd 
12001 iBlitz fc Rosolowskvll200a iBigiel et alJl2008l ). The 
observed trend of the atomic-to-molecular transition col- 
umn density to decrease with increasing metallicity is 
also quantitatively reproduced in our calculations. 

Figure [5] does not illustrate the relationship between 
the atomic and molecular gas inside mostly molecular 



gas, when /h 2 ~ 1- I n order to test our model in that 
regime, we show in Figure [7] a complement of Figure [SJ 
the atomic hydrogen fraction as a function of the total 
gas column density. While the measurements of atomic 
hydroge n fraction in the mostly molecular gas are n ot nu- 
merous dGoldsmith fe Lil 12003: ILerov et all 12001 . they 
provide a useful constraint for our model (see also the 
Appendix) in the regime where the gas is mostly molec- 
ular (/hi < 1). 

Figure [8] shows the dependence of the transition from 
atomic to molecular phase on the adopted clumping fac- 
tor. Although the dependence of the transition on the 
clumping factor 6 is weaker than on metallicity, it is still 
appreciable. Increasing the clumping factor steepens the 
transition and shifts it to lower number densities. Com- 
parison with observations in Figure [6] indicates that the 
clumping factor of C p — 10 is a good fiducial value. This 
value is also within the range of clumping factors ex- 
pected for log-normal PDFs in the theoretical models of 
molecular clouds (see § 12. 2[) . 

In Figure [9] we examine the dependence of the tran- 
sition from atomic to molecular gas on the UV flux. 
The right panel (simulation D) represents a UV flux 
10 times greater than the one on the left panel (sim- 
ulation A); that value is at the upper range of obser- 
vation al estimates of th e UV field in gamma-ray burst 
hosts ()Chen et al.| [2008). Although some differences can 
be seen, they are much smaller than in the previous fig- 
ures. The metallicity of the gas, rather than the amount 
of dissociating radiation, is the primary factor determin- 
ing the density of atomic to molecular transition. This 
result may appear counter- intuitive at first, but it is not 
surprising. Since the UV flux is high, in the optically thin 
regime the equilibrium fraction of molecular hydrogen is 
of the order of 10~ 8 . Thus, gas becomes mostly molec- 
ular only when shielding optical depth is of the order of 

— ln(10 -8 ) « 16. At this high column density, a factor 
of 10 increase in the UV flux (and, thus, a factor of 10 
decrease in the equilibrium H2 fraction in the optically 
thin regime) can be fully compensated by the increase 
of the optical depth (and, thus, the column density) to 

— ln(10 -9 ) « 18, which can be achieved, for example, by 

6 And hence on the unresolved structure of the gas within our 
resolution elements. 
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Fig. 6. — Molecular hydrogen fraction in level 9 cells, plotted against the total neutral gas column density. The grayscale corresponds 
to the density of binned points on the plane. Left panel: simulation A (solar metallicity); middle panel: simulation Bl (metallicity 0.3 
solar); right panel: simulation B 2 (metallicity . 1 sola r). Different symbols sho ws observational data, as follows: open ci rcles: FUSE 
SMC measurements (compiled by Browning ct al. (2003) and Gillmon et al. (2006)); filled diamonds: SMC measurements from Ler ov et al.l 
(2007); open diamonds: FUSE LMC measurements; o pen triangles: FUSE MW Halo measurements; open squares: FUSE MW disk 
measurements; filled squares: MW measurements from [Wolfirc ct al. (200^). The black/yellow striped line in the left panel shows the 
location of points with Shi = IOMq/pc 2 . 

(SF2) and simulation E with a constant star formation 
timescale (SF1). The figure shows that there is no dis- 
cernible differences between these two recipes (we have 
also checked that there is no difference when SF3 recipe 
is adopted). The details of the star formation recipe 
adopted thus do not affect the transition of atomic to 
molecular phase, at least within the range of recipes and 
parameters we considered. 

3.2. The Kennicutt- Schmidt relation 

In this section we examine the large-scale relation 
between star formation and gas surface densities (the 
Kennicutt-Schmidt relation) in the simulated galaxies 
with the molecular gas model and star formation pre- 
scriptions discussed in the previous sections. Figures [TT] 
- Q3] show this relation for different choices of the model 
parameters and star formation prescriptions. Similarly 
to observational estimates, the star formation rate in this 
figure was averaged over 30 Myr, the time interval com- 
parable to the characteristic life time of molecular clouds 
and massive stars. Each point in the figures corresponds 
to a level 3 cell in the simulation at z — 4; the physical 
cell size is 3.3 kpc and thus both the star formation and 
gas surface densities have been averaged on this scale. 
In addition to the total gas surface density, the figures 
also show the dependence of star formation on the sur- 
face density of atomic and molecular gas separately. The 
solid line in each panel corresponds to the best fit to the 
empirical correlat ion estimated for nearby massive and 
starburst galaxies iKennicutti (|1998af) : 
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Fig. 7. — Atomic hydrogen hydrogen fraction in level 9 cells, 
plotted against the total gas number density (a complementary 
figure to Fig. [B). Two panels show solar and 10% solar metallicity 
cases (we have found no observational data for the LMC case, so we 
do not show 30% solar metallicity case). Open and filled triangle for 
the so lar metallicity case show measurements from Goldsmith & Li 
(2005) , while filled ' squares in the 10% solar metallicity panel show 
I Lerov c t al. ( 2007) data. Solid lines in both panels show constraints 
S HI = 10M Q /pc 2 and 100M o /pc 2 respectively. 



only w 15% change in the gas density. 

This is not to say that the UV flux is not at all im- 
portant. Recall that the mean flux within our model 
ISM is already quite high (Figure [2]) . Without this ra- 
diation the ab undance of H2 at smaller den sities would 
be quite high ([Robertson fc Kravtsovl[2008l ) and transi- 
tion to fully molecular phase not as steep as indicated by 
observations (Figure [6]). The role of the interstellar UV 
flux therefore appears to be in controlling the amount 
of diffuse H2 and maintaining the thermodynamic bal- 
ance of the cold neutral medium at lower densities, while 
the transition to the fully molecular phase is controlled 
primarily by metallicity and dumpiness of the gas. 

For completeness, Figure [TU] shows dependence of the 
atomic to molecular phase transition on the adopted star 
formation recipe. In particular, we compare our fiducial 
simulation (A) with a star formation timescale which 
scales with the inverse square root of the gas density 



Ssfr = 2.5x10" 



J gas 



1.4 



M Q yr" 1 kpc" 2 . (8) 



,1M / P c 2 , 

Before we discuss these figures further, it is worth ex- 
plaining certain features of these plots which may look 
peculiar to the reader. When one compares the KS rela- 
tions with different H2 and SF model parameters, Ssfr 
sometimes does not change, while Ehi and Eh 2 change 
dramatically. This seemingly strange behavior is partic- 
ularly apparent at high surface densities. The reason is 
that star formation rate is averaged over 30 Myr, while 
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Fig. 8. — Molecular hydrogen fraction in level 9 cells, plotted against the total gas number density. Panels from the left to right correspond 
to clumping factor C p = 1 (simulation CI), C' p = 10 (simulation A), and C p = 100 (simulation C2) respectively. The grayscale corresponds 
to the density of binned points on the plane. 
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Fig. 9. — Molecular hydrogen fraction in level 9 cells, plotted against the volume number density of the total gas. The right panel 
corresponds to a simulation with 10 times the UV flux (simulation D) of that on the left (simulation A). The grayscale corresponds to the 
density of binned points on the plane. 

the fraction of H2 and H I can change on a shorter time 
scale. Thus, if a region of high surface density becomes 
fully molecular and forms a population of stars, which 
then dissociate much of the remaining H2, Ssfr will re- 
flect the fully molecular E gas and not the current instan- 
taneous Eh 2 - The latter reflects the speed with which 
the molecular gas is able to regenerate after being dis- 
sociated by young stars. This speed is controlled by the 
metallicity and dumpiness of the medium and thus the 
instantaneous £h 2 will reflect these dependencies, while 
Esfr will not. 

Figure QT] shows that the KS relation is quite sensi- 
tive to the metallicity of the gas. The left panel corre- 
sponds to our fiducial simulation (A) with solar metallic- 
ity and the right panel to simulation B2 with metallicity 
0.1 solar. While the star formation rate is tightly corre- 
lated with the total and H2 surface densities, the corre- 



lation with the HI density exhibits considerable scatter. 
In addition, the HI surface density saturates at a rela- 
tively low surface density, which depends on metallicity 
(at « 30M Q /pc 2 for Z = Z Q and « 100M Q /pc 2 for 
Z = O.IZq). S uch saturation is consistent with observa- 
tions (see, e.g-lKennicutt et alJl2007l : iBigiel et~alll2008t 
iLerov et al.ll200ct ). The surface density at which H I satu- 
rates and its metallicity dependence are of course related 
to the volume density of the transition from the atomic 
to molecular phase (see Fig. EJ. Given that the latter 
depends on the choice of the model parameters, the sat- 
uration in the KS relation can be used as an additional 
observational constraint. 

The relationship between the star formation rate and 
the molecular gas surface density ("molecular KS rela- 
tion" ) is less affected by the gas metallicity. Indeed, at a 
given surface density, the star formation rate is about 3 
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Fig. 10. — Molecular hydrogen fraction in level 9 cells, plotted against the total gas number density. Left panel: simulation A (star 
formation recipe SF2); right panel: simulation E (star formation recipe SF1). The grayscale corresponds to the density of binned points 
on the plane. 




Fig. 11. — Sensitivity of the Esfr — S gas relation to metallicity. Quantities are averaged over level 3 cells (corresponding, at z = 4, to a 
proper scale of 3.2 kpc). Star formation rates are averaged over 30 Myrs. Each point corresponds to a different level 3 cell in the simulation 
at 2 = 4. Left panel: simulation A (solar metallicity); right pa nel: simulat i on B2 (metallicity 0.1 solar). Red triangles: molecular gas; blue 
squares: atomic gas; black stars: total gas (no He). Solid line: Kennicutt (1998a) law (slope equal to 1.4, see Eq. |8j. 



times higher in a low metallicity gas. This simply reflects 
the fact that, with our assumptions of a linear relation 
between the dust-to-gas ratio and gas metallicity gas at 
0.1 solar metallicity forms stars at typical densities that 
are 10 times higher than in a solar metallicity gas, result- 
ing in about a factor of 3 higher specific star formation 
rate for our fiducial star formation recipe SF2; respec- 
tively, there would be almost no change in the molecular 
KS relation for the SF1 recipe. 
Figure [12] compares the KS relations for the calcula- 



tions with different clumping factors, C p . Our fiducial 
simulation A (C p — 10) is shown in the middle panel, 
while a simulation with smaller and greater clumping fac- 
tors (simulation CI with C p = 1 and simulation C2 with 
C p = 100) are shown in the left and right panels respec- 
tively. Qualitatively, the results are similar to those in 
Figure[TlJ However, the figure shows significant sensitiv- 
ity of the H I saturation surface density to the clumping 
factor. Again, this is related to the corresponding de- 
pendence of the local volume density of the transition 
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Fig. 12. — Sensitivity of the SgpR — S gas relation to clumping factor. Points and lines as in Fig. 1111 Left panel: C' p = 1 (simulation CI); 
middle panel: C p = 10 (simulation A); right panel: C' p = 100 (simulation C2). 
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Fig. 13. — Sensitivity of the SgpR — S gaa relation to the UV flux. Points and lines as in Fig. 1111 The right panel corresponds to a 
simulation with 10 times the UV flux (simulation D) than that on the left (simulation A). 



shown in Figure [8j The H I saturation surface density 
observed in nearby galaxies can place constraints on the 
value of the clumping factors. For example, the gas in 
M51 has metallicity of Z ss 2Zq and its KS relation ex- 
hibits saturation at Ehi ~ 30M Q /pc 2 . In a larger sam- 
ple of both massive spirals a nd dwarf gal a xies ( with typ- 
ically Z < Z Q ) reported by iLerov et all (|2008h . the HI 
saturation occurs at Shi ~ 10 — 2OM /pc 2 . Comparing 
these numbers to the KS relations in Figure [T2l indicates 
that the clumping factor required to reproduce these ob- 
servations (especially at lower than solar metallicities) is 
quite high, C p > 10. This is, in fact, consistent with ob- 
servational estimates of the clumping within molecular 
clouds which indicate that typical ratios of the local to 
the mean density of a clo ud are ~ 30 — 100 (see §3.1.2 in 
iMcKee fc Ostriker1l2007L and references therein). 

As Figure[8jdemonstrates, the column density at which 



the H I to H2 transition occurs decreases with increasing 
clumping factor; as a result, the surface density at which 
HI saturates decreases with increasing clumping factor. 
Consequently, the star formation rate scaling with the 
H2 gas approaches the scaling with the total gas as the 
clumping factor increases, and become very similar for 
the highest clumping factor we have used (C p = 100). 

The sensitivity of the KS relation to the UV flux is 
shown in Figure 1131 The right panel corresponds to a 
simulation with 10 times the UV flux (simulation D) than 
that that of our fiducial run on the left (simulation A). 
No appreciable changes can be seen between the two pan- 
els. The dependence of the global star formation on the 
UV flux is therefore much weaker than on the metallicity 
and dumpiness. As we have noted above, this does not 
mean that the UV flux is not important at all. With the 
UV flux close to zero the abundance of H2 at smaller den- 



11 




io" 2 io _1 i 10 1 io 2 io 3 io~ 2 io _1 i io 1 io 2 io 3 

2 gas (M pc" 2 ) E gas (M pc" 2 ) 



Fig. 14. — Sensitivity of the Esfr — E gas relation to the star formation recipe. Points and lines as in Fig. 1111 Left panel: star formation 
recipe SF2 (simulation A); right panel: star formation recipe SF1 (simulation E). To illustrate the difference in the slope better, we also 
show a KS relation with the slope of 1 for S > 30Mq/pc 2 . with a dotted line. 



sities would be quite high (|Robertson fe Kravtsovll2008f ) 
and transition to fully molecular phase not as steep as 
indicated by observations (Figure E]). The role of the in- 
terstellar UV flux therefore appears to be in controlling 
the amount of diffuse H2 and maintaining the thermody- 
namic balance of the cold neutral medium at lower den- 
sities. Beyond certain level, however, the dependence of 
the results on the UV flux saturates. 

Finally, Figure [14] shows the dependence of the scal- 
ing on the local star formation recipe used. The left 
panel corresponds to our fiducial simulation (A), with 
star formation recipe SF2 (timescale scaling with den- 
sity), while the right panel corresponds to simulation E 
with star formation recipe SF1 (constant timescale). The 
main difference is the slope of the Ssfr — Sjj 2 relation. 
For the constant time scale recipe SF1 the slope is shal- 
lower and is close to unity. Note, however, that the slope 
of the Ssfr — S gas relation does not change as dramat- 
ically and in fact is fairly c onsistent with the best fit 
relation of iKennicuttJ |l998b). This is because this slope 
is controlled both by the dependence of local star forma- 
tion efficiency on density and on the de pendence of mas s 
fraction of star forming regions on S gas (|Kravtsovll2003h . 

4. DISCUSSION AND CONCLUSIONS 

Specifics of star formation modeling in cosmologi- 
cal simulations of galaxy formation have critical im- 
pact on the properties of resulting galaxies, such as 
their star formation histories (and hence luminosity 
and colors) and morphology. Although significant suc- 
cesses have been achieved using s uch high- resolution 
simulations (c.f., Ma ver et al.l [2008. for a comprehen- 
sive review), many challenges remain. In particular, 
the two likely related challenges arc the low efficiency 
of star formation in low mass systems and the preva- 
lence of thin disks among galaxies. So far, the an- 
swer for these challenges was to devise efficient schemes 
of stellar energy feedback, which is capable of driving 



outflows fe.g.. [Springel k, Hernquistll2003t IStinson et al.1 
200G; C everino &: Klvpinll2007f l. However, at least part 
of the solution may be due to the inh erently low efficiency 
of ga s conversion into stars (e.g., iKrumholz fc McKeel 
2005). Indeed, both nearby and high-redshift galaxies 
provide a variety of clues suggesting that gas conversion 
into stars in low-mass, low-metallicity galaxies is very 
inefficient. To explore the effects of such inefficiency on 
galaxy evolution, a more realistic way of identifying and 
treating star forming regions is needed. In particular, 
the global efficiency with which a galaxy converts gas 
into stars depends on its ability to convert a significant 
fraction of its gas mass into fully molecular form, within 
which conditions for star formation can be realized. The 
latter depends on the small-scale gas properties and local 
interstellar radiation field. These dependencies are not 
captured in the standard recipes of star formation. 

In this paper we have presented a model for molecu- 
lar hydrogen formation for cosmological galaxy formation 
simulations designed to follow the transition from the 
atomic to molecular phase on the scales of tens of parsecs. 
The model is applicable in simulations, in which individ- 
ual star forming regions - the giant molecular complexes 
- can be identified and their mean internal density esti- 
mated reliably, even if internal structure is not resolved. 
We present a number of tests of the model and illustrate 
its effects on the global correlation between star forma- 
tion rate and gas surface density. 

The model shows that the transition from atomic to 
fully molecular phase depends primarily on the metallic- 
ity (which we assume is directly related to the dust abun- 
dance) and dumpiness of the interstellar medium. The 
dumpiness simply boosts the formation rate of molec- 
ular hydrogen, while dust both serves as a catalyst of 
H2 formation and as additional shielding from dissoci- 
ating UV radiation. The upshot is that it is difficult 
to form fully-shielded giant molecular clouds, while gas 
metallicity is low. However, once the gas is enriched 
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(say to Z ~ 0.01 — O.IZq), the subsequent star forma- 
tion and enrichment can proceed at a much faster rate. 
This may keep star formation efficiency in the low-mass, 
low-metallicity progenitors of galaxies very low for a cer- 
tain period of time. One can think of this as a "feed- 
back" mechanism because it reduces star formation in 
low metallicity galaxies, but also speeds up star forma- 
tion once conditions for sustained conversion of a large 
fraction of gas into fully molecular form are realized. 

Molecular fractions in the low-metallicity environ- 
ments of the nearby dwarf and low surface brightness 
galaxies are ind eed very small ([Matthews et al.1 120051 : 
iDas et al\ l2006f ) . Recent observational evidence indi- 
cates that the fraction of star forming molecular gas is 
also small in high-redshift g alaxies (e.g.. lTumlinson et al.l 
120071 1 Wolfe fc Chenl 12009) . Dependence of the density 
at which transition from atomic to fully molecular phase 
occurs on metallicity is indeed observed both in direct 
measurements for nearby galaxies (see Fig. [6J and in 
the maximum HI column den sities of DLA s of different 
metallicity at high redshifts (|Schavd I2001D and in the 
surface density of gas at which surface density of HI sat- 
urates in nearby galaxies (jKrumholz et al.ll2008b| ). 

The UV radiation field is also important in main- 
taining the cold neutral m edium in atomic form 
([Robertson fc Kravtsovll2008l ). Our results show, how- 
ever, that beyond certain flux level the effect of UV ra- 
diation on the ability of ISM to form molecular clouds 
saturates. 

We show that the global Kennicutt-Schmidt relation 
between gas and star formation surface densities is also 
affected by the local processes controlling conversion of 
atomic gas into molecu lar. Our results are qualitativ ely 
consistent with those of lRobertson fc Kravtsovi (|2008l ) in 
that star formation rate dependence on the total gas den- 
sity can vary depending on the local conditions in the 
interstellar medium. These results are also consistent 
with existing detailed studies of the Kennicutt-Schmidt 
relation in nearby gala xies with lower metallicities and 
surfa ce densities (e.g., iBoissier et al.l [20031 : iHeyer et al.1 
12004 Wyder et al. 2008, in preparation). 

The most recent comprehensive st udy of the KS rela - 
tion in the THINGS galaxy sample bv lBigiel et all (j2008l ) 
shows the KS relation qualitatively similar with our re- 
sults presented in Figures [TT] - [111 Esfr dependence 
on E gas is steep at low gas surface densities, but be- 
comes shallow at larger E gas values. Interestingly, they 
find evidence that the slope n of the KS relation at 
E g as ~ 10 - 100 M is n « 1 ± 0.2, shallower than the 
canonical value of n ss 1.4. The relation at higher sur- 
fa ce densities in star burst galaxies steepens (see Fig. 15 
of lBigiel etali r2008). In the context of our calculations, 
these results can be interpreted as the systematic change 
of the internal density of molecular clouds. In our star 
formation recipes the slope of the large-scale relation de- 



pends on the assumed dependence of star formation time 
scale, Tgf , on the local internal density on the scale of star 
forming regions, as shown in Figure 1141 The steepen- 
ing of the observed KS relation at high surface densities 
can thus be interpreted as due to systematic increase of 
the internal density of molecular clouds at high surface 
densities, as can be expected in the high-pressure envi- 
ronm ents of starburst galaxies (e.g. iDownes fc Solomon! 
Il998h . 

Another prediction of our calculations is that the tran- 
sition at low gas surface densities from a steep to a 
shallow relation should depend on metallicity. It would 
thus be interesting to also explore such metallicity de- 
pendence in observations. As a first step, one can con- 
sider comparison of the KS relation in lower-metallicity 
dwarf galaxies to the low-Ey as regions in the ou tskirts of 
higher-metallicity spirals by Bi giel et al.l (|2008l see their 
Fig. 12). Although the KS relation is steep in both of 
these regimes, the relation for the outskirts of spirals is 
shifted somewhat to lower E gas , as could be expected 
from the metallicity dependence. 

The above comparisons illustrate that a star forma- 
tion model of the kind prese nted in this paper (see also 
iRobertson fc Kravtsovll2008f) can be very useful in inter- 
preting the detailed features of the observed KS relation. 
Conversely, one can use the detailed observations of the 
KS relation in different regimes (different E gas , different 
metallicities, etc.) to constrain and tune model param- 
eters. This, in turn, may significantly improve fidelity 
of star formation modeling in simulations. We plan to 
present such comparisons as well as applications of our 
model in fully self-consistent cosmological simulations of 
galaxy formation in a forthcoming work. 
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APPENDIX 
Numerical Considerations 

Ideally, in a simulation with infinite spatial resolution, our model should work as designed. In practice, however, 
the spatial resolution of a simulation is finite. In addition, numerical solutions of partial differential equations always 
contain truncation errors. In our case, these errors lead to a small, but non-negligible amount of numerical heating 
and advection that biases the H2 formation model. 

An example of such diffusion is shown in Figure 1151 where we compare an ideal, single zone (i.e. without any 
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Fig. 15. — Evolution of the gas temperature (top) and HI fraction (bottom) for a single, fully molecular cell in our fiducial run. The 
black dashed line shows a test calculation when the cell is evolved in isolation at constant gas density and constant radiation field (without 
any hydrodynamic effects). The light gray solid line is the evolution of the same cell in a full cosmological simulation - notice the numerical 
heating and advection, appearing as jumps in temperature and neutral hydrogen fraction at regular intervals corresponding to individual 
hydrodynamic time steps. The solid black line is for the same test with the numerical correction l|3) l. 



hydrodynamic effects) evolution of a typical cell inside molecular clouds (i.e. with molecular fraction above 99% - 
shown as black dashed lines) with the evolution of this cell in a real cosmological simulation. Because our resolution 
is finite, gas motions across individual cells are not fully resolved; in a Riemann code the truncation errors from such 
incomplete resolution lead to cells being heated at each time step. This numerical heating is not too serious, since 
the gas cooling time is very short; however, a more serious artifact is the numerical advection of atomic gas inside 
molecular clouds. The gray line on the bottom panel shows the evolution of the atomic hydrogen fraction in that cell 
- at each hydrodynamic time step there is a small (about 10%) increase in the atomic hydrogen fraction. This amount 
is not large, but it happens at every time step, and so, effectively, there is a numerical flux of atomic hydrogen into 
molecular clouds, with the corresponding opposite flux of molecular hydrogen out of molecular gas. 

We emphasize that this is a numerical artifact: the transition between molecular and atomic phases is so sharp that 
few numerical schemes are able to treat it perfectly. As a numerical artifact, this effect is actually rather small - 10% 
jump in /hi ~ 0.01 corresponds to the relative numerical diffusion of only about 10~ 3 per time step. This effect is 
only important because the atomic hydrogen fraction inside molecular clouds is of the same order, about 10 -3 , or even 
less. 

In order to mitigate this numerical advection, we adopt the following, entirely empirical approach. The numerical 
advection can be thought of as an additive increase in the atomic hydrogen fraction per time step, 

HI — A HI + AA HI- 
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If we multiply the "numerical" value Xg j m at the end of each hydrodynamic time step by the correction factor 

JNum — ^1 "I" j^Truc 

then we would remove the numerical advection completely. Unfortunately, we do not have any way to measure both 
Xjjj uc and AXhi in a self-consistent way (for example, keeping track of -Xhi without hydrodynamic advection would 
violate Galilean invariance). Because the numerical advection per time step is not large, we can use X^j m (the value 
we get in the simulation) instead of ^hi uc m tne expression for fy um , but the factor AXhi needs to be deduced 
heuristically. 
We choose to define AXh i as 

AX H I = aa- C ell-r— X H , (1) 



A 



x 



where a is a numerical coefficient of the order of unity, <7 ce n is the gas velocity dispersion at the cell scale, computed 
in each cell k) from its six neighbors, 

0-cell(M! k ) = g (Vi+l,j,k ~ V l}j , k f + {Vi-l,j,k ~ Vi,],k) 2 + (#i,j+l,k - V^ hk f + ... , (2) 

and At and Ax are the numerical time step and the cell size (which is, in fact, different for different cells on an 
adaptively refined mesh in our simulations). The factor Xjj encapsulates our assumption that the gas is mostly atomic 
and neutral just outside molecular clouds. 

Equation ([I]) has the correct properties for an expression describing numerical advection. First, it vanishes in the 
limit of At — > 0. We have actually verified, by running test simulations with reduced time steps, that the numerical 
advection effect shown in Figure [15] scales linearly with the time step At. Second, equation [T] is Galilean invariant 
and depends on the local variation in the gas velocity, i.e. on the quantity that determines hydrodynamic advection 
(beyond trivial uniform translation). Third, AXhi does not actually diverge in the limit Ax — > 0, because for a fully 
resolved flow 

1 / dv l dv l \ 



dxi dxi 

Thus, our heuristic numerical advection correction factor becomes 

/ At XnY 1 

/Nu m ^l + aa cell — — j . (3) 

This form for the numerical correction is also rather insensitive to the specific choice of the coefficient a: varying a 
from 0.5 to 2 changes atomic hydrogen fractions inside molecular clouds and star formation rates by less than their 
natural scatter. 

The correction factor ([3]), being heuristic, does not completely correct for numerical advection on a cell- by-cell 
basis; it is valid only in a statistical sense. As Figure [T5] shows, evolution of the abundance in a test cell with the 
correction factor included does not match a single zone calculation (although, a single zone calculation excludes not 
only unphysical numerical advection, but physical real advection as well, and should not be taken as a correct solution). 
The only justification for the factor ([3|) is Figure [TBI which sho ws a comparison of the a tomic hydrogen fractions that 
we find in our simulations with the observational points from iGoldsmith fc Lil (|2005l ). Without the correction {3} 
the observational points cannot be reproduced, while with the correction we get a good match to the observational 
constraints. 
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